Machine Learning - Homework 2


In this exercise, we deal with crimes data obtained by the Isareli police.

Using various machine learning algorithms, we analyze the data to find patterns and predict future crime rates.


0 Imports

Note: we used plotly in this project. Thus, make sure to install it first.

In [1]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.impute import KNNImputer
from sklearn.preprocessing import StandardScaler

# question 3
from sklearn.cluster import KMeans
from sklearn.mixture import GaussianMixture
from sklearn.decomposition import PCA
import plotly as py
import plotly.graph_objs as go
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot

# question 4
from sklearn.model_selection import GridSearchCV
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error, confusion_matrix, accuracy_score

# question 5
from sklearn.ensemble import AdaBoostRegressor
from matplotlib.pyplot import *

# Set random seed to ensure reproducible runs
RSEED = 100

1 Data Cleaning and Wrangling

The data include three files.

  1. Crimes - the number of crimes of different categories for each city, throughout the years 2014-2019.
  2. Cities - demographic data of the cities.
  3. Mapping - maps the different categories of the cities file.

1.1 Crimes Dataset

In [2]:
crimes = pd.read_excel('crimes.xlsx')
crimes
Out[2]:
סמל יישוב יישוב מחושב תאור קבוצה סטטיסטית שנת הודעה Total 2019 2018 2017 2016 2015 2014
0 NaN Total NaN NaN 1973220 301142 320713 329265 328681 339804 353615
1 472.0 אבו גוש Total NaN 1997 338 284 310 340 323 402
2 472.0 אבו גוש - NaN 4 2 2 - - - -
3 472.0 אבו גוש עבירות בטחון NaN 55 7 9 8 6 15 10
4 472.0 אבו גוש עבירות כלכליות NaN 3 1 - - - 2 -
... ... ... ... ... ... ... ... ... ... ... ...
3678 1054.0 תל שבע עבירות נגד גוף NaN 1107 179 188 202 189 171 178
3679 1054.0 תל שבע עבירות סדר ציבורי NaN 2112 427 444 344 316 262 319
3680 1054.0 תל שבע עבירות רשוי NaN 41 6 10 9 8 1 7
3681 1054.0 תל שבע עבירות תנועה NaN 67 7 13 16 12 9 10
3682 1054.0 תל שבע שאר עבירות NaN 7 - 3 2 1 - 1

3683 rows × 11 columns

Two column names are changed for convenience.

In [3]:
names_dict = {'יישוב מחושב': 'שם יישוב', 'תאור קבוצה סטטיסטית': 'עבירה'}
crimes.rename(columns=names_dict, inplace=True)
crimes.head()
Out[3]:
סמל יישוב שם יישוב עבירה שנת הודעה Total 2019 2018 2017 2016 2015 2014
0 NaN Total NaN NaN 1973220 301142 320713 329265 328681 339804 353615
1 472.0 אבו גוש Total NaN 1997 338 284 310 340 323 402
2 472.0 אבו גוש - NaN 4 2 2 - - - -
3 472.0 אבו גוש עבירות בטחון NaN 55 7 9 8 6 15 10
4 472.0 אבו גוש עבירות כלכליות NaN 3 1 - - - 2 -

The column שנת הודעה includes only one unique value and thus will be dropped.

In [4]:
len(crimes['שנת הודעה'].unique())
Out[4]:
1
In [5]:
crimes.drop(['שנת הודעה'], axis=1, inplace=True)

In addition, all rows and columns that are a sum of others, will be dropped because of innacurate values.

In [6]:
crimes.drop(['Total'], axis=1, inplace=True)
crimes = crimes.loc[crimes['שם יישוב'] != 'Total']
crimes = crimes.loc[crimes['עבירה'] != 'Total']
In [7]:
crimes.head()
Out[7]:
סמל יישוב שם יישוב עבירה 2019 2018 2017 2016 2015 2014
2 472.0 אבו גוש - 2 2 - - - -
3 472.0 אבו גוש עבירות בטחון 7 9 8 6 15 10
4 472.0 אבו גוש עבירות כלכליות 1 - - - 2 -
5 472.0 אבו גוש עבירות כלפי המוסר 14 5 11 6 8 16
6 472.0 אבו גוש עבירות כלפי הרכוש 93 91 94 95 87 115

Next, all missing values that represent the number of crimes (i.e. columns 2014-2019), are replaced with 0.

In [8]:
years = [2014, 2015, 2016, 2017, 2018, 2019]
crimes[years] = crimes[years].replace('-', 0)
In [9]:
crimes.head()
Out[9]:
סמל יישוב שם יישוב עבירה 2019 2018 2017 2016 2015 2014
2 472.0 אבו גוש - 2 2 0 0 0 0
3 472.0 אבו גוש עבירות בטחון 7 9 8 6 15 10
4 472.0 אבו גוש עבירות כלכליות 1 0 0 0 2 0
5 472.0 אבו גוש עבירות כלפי המוסר 14 5 11 6 8 16
6 472.0 אבו גוש עבירות כלפי הרכוש 93 91 94 95 87 115

In some cities, there are crimes of category לא ידוע or a missing value (-). These will all be merged into the category שאר העבירות as they are not specific to a known category.

In [10]:
names = crimes['שם יישוב'].unique()

for name in names:
    
    # get rows of specific city
    x = crimes.loc[crimes['שם יישוב'] == name]
    
    # get values of crime category '-'
    x1 = x.loc[(x['עבירה'] == '-'), years]
    # get values of crime category 'rest'
    x2 = x.loc[x['עבירה'] == 'שאר העבירות', years]
    # get values of crime category 'unknown'
    x3 = x.loc[x['עבירה'] == 'לא ידוע', years]
    
    # sum of values
    s = x1.values.astype(int) + x2.values.astype(int) + x3.values.astype(int)
    # append sum to df
    crimes.loc[x2.index.values, years] = s

# drop values of crime category '-' and unknown
crimes = crimes.loc[(crimes['עבירה'] != '-') & (crimes['עבירה'] != 'לא ידוע')]
In [11]:
crimes.head(12)
Out[11]:
סמל יישוב שם יישוב עבירה 2019 2018 2017 2016 2015 2014
3 472.0 אבו גוש עבירות בטחון 7.0 9.0 8.0 6.0 15.0 10.0
4 472.0 אבו גוש עבירות כלכליות 1.0 0.0 0.0 0.0 2.0 0.0
5 472.0 אבו גוש עבירות כלפי המוסר 14.0 5.0 11.0 6.0 8.0 16.0
6 472.0 אבו גוש עבירות כלפי הרכוש 93.0 91.0 94.0 95.0 87.0 115.0
7 472.0 אבו גוש עבירות מין 5.0 9.0 4.0 7.0 5.0 3.0
8 472.0 אבו גוש עבירות מרמה 5.0 9.0 7.0 10.0 14.0 16.0
9 472.0 אבו גוש עבירות נגד אדם 0.0 0.0 2.0 3.0 1.0 2.0
10 472.0 אבו גוש עבירות נגד גוף 113.0 105.0 123.0 138.0 116.0 119.0
11 472.0 אבו גוש עבירות סדר ציבורי 202.0 167.0 175.0 191.0 185.0 251.0
12 472.0 אבו גוש עבירות רשוי 0.0 0.0 0.0 0.0 1.0 7.0
13 472.0 אבו גוש עבירות תנועה 1.0 0.0 2.0 4.0 2.0 5.0
14 472.0 אבו גוש שאר עבירות 0.0 0.0 1.0 0.0 0.0 2.0

The crime rates throughout the years are highly correlated. Therefore, we will use the average of years 2014-2018 for our train dataset, and the year 2019 for our test dataset. We've decided to use the average because it's a suffecient statistic and an unbiased estimate for the number of crimes.

In [12]:
crimes[years].corr()
Out[12]:
2014 2015 2016 2017 2018 2019
2014 1.000000 0.997521 0.992773 0.989606 0.982249 0.965137
2015 0.997521 1.000000 0.995270 0.990489 0.982768 0.967562
2016 0.992773 0.995270 1.000000 0.994767 0.989679 0.977627
2017 0.989606 0.990489 0.994767 1.000000 0.996011 0.983887
2018 0.982249 0.982768 0.989679 0.996011 1.000000 0.991187
2019 0.965137 0.967562 0.977627 0.983887 0.991187 1.000000
In [13]:
# add average column
cols = crimes.loc[:, years[0:5]]
crimes['ממוצע'] = cols.mean(axis=1)
In [14]:
crimes.head(12)
Out[14]:
סמל יישוב שם יישוב עבירה 2019 2018 2017 2016 2015 2014 ממוצע
3 472.0 אבו גוש עבירות בטחון 7.0 9.0 8.0 6.0 15.0 10.0 9.6
4 472.0 אבו גוש עבירות כלכליות 1.0 0.0 0.0 0.0 2.0 0.0 0.4
5 472.0 אבו גוש עבירות כלפי המוסר 14.0 5.0 11.0 6.0 8.0 16.0 9.2
6 472.0 אבו גוש עבירות כלפי הרכוש 93.0 91.0 94.0 95.0 87.0 115.0 96.4
7 472.0 אבו גוש עבירות מין 5.0 9.0 4.0 7.0 5.0 3.0 5.6
8 472.0 אבו גוש עבירות מרמה 5.0 9.0 7.0 10.0 14.0 16.0 11.2
9 472.0 אבו גוש עבירות נגד אדם 0.0 0.0 2.0 3.0 1.0 2.0 1.6
10 472.0 אבו גוש עבירות נגד גוף 113.0 105.0 123.0 138.0 116.0 119.0 120.2
11 472.0 אבו גוש עבירות סדר ציבורי 202.0 167.0 175.0 191.0 185.0 251.0 193.8
12 472.0 אבו גוש עבירות רשוי 0.0 0.0 0.0 0.0 1.0 7.0 1.6
13 472.0 אבו גוש עבירות תנועה 1.0 0.0 2.0 4.0 2.0 5.0 2.6
14 472.0 אבו גוש שאר עבירות 0.0 0.0 1.0 0.0 0.0 2.0 0.6
In [15]:
# create train dataset
crimes_train = crimes[['סמל יישוב', 'שם יישוב', 'עבירה', 'ממוצע']]
crimes_train.head()
Out[15]:
סמל יישוב שם יישוב עבירה ממוצע
3 472.0 אבו גוש עבירות בטחון 9.6
4 472.0 אבו גוש עבירות כלכליות 0.4
5 472.0 אבו גוש עבירות כלפי המוסר 9.2
6 472.0 אבו גוש עבירות כלפי הרכוש 96.4
7 472.0 אבו גוש עבירות מין 5.6
In [16]:
# create test dataset
crimes_test = crimes[['סמל יישוב', 'שם יישוב', 'עבירה', 2019]]
# rename column 2019
crimes_test = crimes_test.rename({2019: 'מספר עבירות'}, axis='columns')
crimes_test.head()
Out[16]:
סמל יישוב שם יישוב עבירה מספר עבירות
3 472.0 אבו גוש עבירות בטחון 7.0
4 472.0 אבו גוש עבירות כלכליות 1.0
5 472.0 אבו גוש עבירות כלפי המוסר 14.0
6 472.0 אבו גוש עבירות כלפי הרכוש 93.0
7 472.0 אבו גוש עבירות מין 5.0

For easier implementation, we reshape the dataset to a wide format where each row represents one city and the columns are the different types of crimes.

In [17]:
crimes_train_wide = crimes_train.pivot_table(index=['שם יישוב', 'סמל יישוב'], columns=['עבירה'], values=['ממוצע'])
crimes_train_wide.head()
Out[17]:
ממוצע
עבירה סעיפי הגדרה עבירות בטחון עבירות כלכליות עבירות כלפי המוסר עבירות כלפי הרכוש עבירות מין עבירות מנהליות עבירות מרמה עבירות נגד אדם עבירות נגד גוף עבירות סדר ציבורי עבירות רשוי עבירות תנועה שאר עבירות
שם יישוב סמל יישוב
אבו גוש 472.0 NaN 9.6 0.4 9.2 96.4 5.6 NaN 11.2 1.6 120.2 193.8 1.6 2.6 0.6
אבו סנאן 473.0 NaN 12.0 0.2 38.0 77.0 1.2 NaN 5.8 0.2 71.0 133.4 2.0 2.0 0.2
אבן יהודה 182.0 1.0 0.4 0.2 14.8 260.0 7.4 0.4 12.4 0.4 49.2 104.6 0.4 1.2 0.2
אום אל פחם 2710.0 1.2 118.4 4.2 112.0 433.4 7.0 0.4 47.0 7.0 240.6 655.6 5.6 11.6 0.6
אופקים 31.0 1.6 9.0 1.6 119.0 349.0 26.0 0.2 33.4 2.2 223.6 439.2 11.4 7.2 1.4
In [18]:
# replace missing values
crimes_train_wide = crimes_train_wide.replace(np.nan, 0)
In [19]:
# set columns and index
crimes_train_wide.columns = crimes_train_wide.columns.droplevel()
crimes_train_wide.reset_index(inplace=True)
In [20]:
crimes_train_wide.head()
Out[20]:
עבירה שם יישוב סמל יישוב סעיפי הגדרה עבירות בטחון עבירות כלכליות עבירות כלפי המוסר עבירות כלפי הרכוש עבירות מין עבירות מנהליות עבירות מרמה עבירות נגד אדם עבירות נגד גוף עבירות סדר ציבורי עבירות רשוי עבירות תנועה שאר עבירות
0 אבו גוש 472.0 0.0 9.6 0.4 9.2 96.4 5.6 0.0 11.2 1.6 120.2 193.8 1.6 2.6 0.6
1 אבו סנאן 473.0 0.0 12.0 0.2 38.0 77.0 1.2 0.0 5.8 0.2 71.0 133.4 2.0 2.0 0.2
2 אבן יהודה 182.0 1.0 0.4 0.2 14.8 260.0 7.4 0.4 12.4 0.4 49.2 104.6 0.4 1.2 0.2
3 אום אל פחם 2710.0 1.2 118.4 4.2 112.0 433.4 7.0 0.4 47.0 7.0 240.6 655.6 5.6 11.6 0.6
4 אופקים 31.0 1.6 9.0 1.6 119.0 349.0 26.0 0.2 33.4 2.2 223.6 439.2 11.4 7.2 1.4
In [21]:
# reshape test dataset
crimes_test_wide = crimes_test.pivot_table(index=['שם יישוב', 'סמל יישוב'], columns=['עבירה'], 
                                      values=['מספר עבירות'], aggfunc='first')
# replace nan values with 0
crimes_test_wide = crimes_test_wide.replace(np.nan, 0)

# set columns and index
crimes_test_wide.columns = crimes_test_wide.columns.droplevel()
crimes_test_wide.reset_index(inplace=True)

1.2 Cities Dataset

In [22]:
cities = pd.read_excel('cities.xlsx')
cities.head()
Out[22]:
שם יישוב סמל יישוב תעתיק מחוז נפה אזור טבעי מעמד מונציפאלי שיוך מטרופוליני דת יישוב סך הכל אוכלוסייה 2018 ... שנת ייסוד צורת יישוב שוטפת השתייכות ארגונית קואורדינטות גובה ועדת תכנון מרחב משטרה שנה שם יישוב באנגלית אשכול רשויות מקומיות
0 אבו ג'ווייעד (שבט) 967 ABU JUWEI'ID 6 62 623.0 NaN NaN 3.0 NaN ... NaN 460 NaN 2.040057e+09 NaN 699.0 15003711.0 2018 Abu Juway'ad NaN
1 אבו גוש 472 ABU GHOSH 1 11 111.0 99.0 444.0 2.0 7543.0 ... NaN 280 NaN 2.105263e+09 598.0 152.0 10002475.0 2018 Abu Ghosh NaN
2 אבו סנאן 473 ABU SINAN 2 24 245.0 99.0 NaN 2.0 13915.0 ... NaN 270 NaN 2.160776e+09 19.0 252.0 10004315.0 2018 Abu Sinan NaN
3 אבו סריחאן (שבט) 935 ABU SUREIHAN 6 62 623.0 NaN NaN 3.0 NaN ... NaN 460 NaN 1.865057e+09 NaN 699.0 10001937.0 2018 Abu Surayhan NaN
4 אבו עבדון (שבט) 958 ABU ABDUN 6 62 623.0 NaN NaN 3.0 NaN ... NaN 460 NaN 1.850058e+09 NaN 699.0 10001937.0 2018 Abu 'Abdun NaN

5 rows × 23 columns

The column שנה includes only one unique value and thus will be dropped.

In [23]:
len(cities['שנה'].unique())
Out[23]:
1
In [24]:
cities.drop(['שנה'], axis=1, inplace=True)

In addition, the columns שם יישוב באנגלית and תעתיק are identical and irrelevant. Therefore, they will be dropped.

In [25]:
cities.drop(['שם יישוב באנגלית', 'תעתיק'], axis=1, inplace=True)

Correlated Features

In [26]:
temp = cities.drop(['שם יישוב'], axis=1)

temp_val=temp.corr().columns.values 
temp_val[7]='סך הכל אוכלוסייה 8102'#to make sure that the reverse will work also to the number 2018
plt.figure(figsize=(15, 12))
sns.heatmap(temp.corr(), cmap='Blues',xticklabels=[i[::-1] for i in temp_val],
            yticklabels=[i[::-1] for i in temp_val])
plt.title('Heatmap of Cities Dataset')
plt.show()

We see that the following columns are highly correlated (a sum of each other), and thus two of them will be dropped.

In [27]:
cols = cities[['סך הכל אוכלוסייה 2018', 'יהודים ואחרים', 'מזה: יהודים', 'ערבים']]
cols.corr()
Out[27]:
סך הכל אוכלוסייה 2018 יהודים ואחרים מזה: יהודים ערבים
סך הכל אוכלוסייה 2018 1.000000 0.965097 0.969768 0.713087
יהודים ואחרים 0.965097 1.000000 0.999038 0.500566
מזה: יהודים 0.969768 0.999038 1.000000 0.523208
ערבים 0.713087 0.500566 0.523208 1.000000
In [28]:
cities.drop(['יהודים ואחרים', 'מזה: יהודים'], axis=1, inplace=True)

In addition, the following columns are also highly correlated, and thus we will only keep the נפה column. Some of the others have missing values and the מחוז column is not as specific as נפה.

In [29]:
cols = cities[['מחוז', 'נפה', 'אזור טבעי', 'ועדת תכנון', 'אשכול רשויות מקומיות']]
cols.corr()
Out[29]:
מחוז נפה אזור טבעי ועדת תכנון אשכול רשויות מקומיות
מחוז 1.000000 0.996007 0.995625 0.992260 0.962393
נפה 0.996007 1.000000 0.999864 0.987192 0.961155
אזור טבעי 0.995625 0.999864 1.000000 0.986150 0.961583
ועדת תכנון 0.992260 0.987192 0.986150 1.000000 0.963143
אשכול רשויות מקומיות 0.962393 0.961155 0.961583 0.963143 1.000000
In [30]:
cities.drop(['מחוז', 'אזור טבעי', 'ועדת תכנון', 'אשכול רשויות מקומיות'], axis=1, inplace=True)

Handling Missing Values

In [31]:
cities.isnull().sum()
Out[31]:
שם יישוב                   0
סמל יישוב                  0
נפה                        0
מעמד מונציפאלי            77
שיוך מטרופוליני          929
דת יישוב                 236
סך הכל אוכלוסייה 2018    264
ערבים                    938
שנת ייסוד                393
צורת יישוב שוטפת           0
השתייכות ארגונית         671
קואורדינטות               32
גובה                     196
מרחב משטרה               164
dtype: int64

For the column השתייכות ארגונית, we assume that the missing values represent 'ללא השתייכות'.

In [32]:
cities['השתייכות ארגונית'] = cities['השתייכות ארגונית'].fillna(19)

A similar approach is used for the column מעמד מונציפאלי where the missing values are assumed to be 'חסר מעמד'.

In [33]:
cities['מעמד מונציפאלי'] = cities['מעמד מונציפאלי'].fillna('חסר מעמד')

For the columns דת יישוב and מרחב משטרה, we added a category of 'אחר' for the missing values.

In [34]:
cities['מרחב משטרה'] = cities['מרחב משטרה'].fillna('אחר')
In [35]:
cities['דת יישוב'] = cities['דת יישוב'].fillna('אחר')

For the column ערבים, if the religion of the city is either 1 - Jewish, or 3 - Bedouin, then the number of arabs is 0.

In [36]:
col = 'ערבים'
cities.loc[cities['דת יישוב'] == 1, col] = cities.loc[cities['דת יישוב'] == 1, col].replace(np.nan, 0)
cities.loc[cities['דת יישוב'] == 3, col] = cities.loc[cities['דת יישוב'] == 3, col].replace(np.nan, 0)

In the column שנת ייסוד there are a few values of 'ותיק' which will be replaced with NaN and later imputed.

In [37]:
cities = cities.replace('ותיק', np.nan)

Next, we decided to drop any rows that have missing values in more than half the columns, and any columns that have missing values in more than half the rows.

In [38]:
cities.shape
Out[38]:
(1482, 14)
In [39]:
cities.dropna(thresh=8, axis=0, inplace=True)
In [40]:
cities.dropna(thresh=800, axis=1, inplace=True)
In [41]:
cities.shape
Out[41]:
(1482, 13)

All the remaining missing values of numeric columns, are imputed using KNN.

In [42]:
imputer = KNNImputer(n_neighbors=5, weights="uniform")

cols = ['סך הכל אוכלוסייה 2018', 'שנת ייסוד', 'קואורדינטות', 'גובה', 'ערבים']

imputed = imputer.fit_transform(cities[cols])
cities[cols] = imputed
In [43]:
cities.isnull().sum()
Out[43]:
שם יישוב                 0
סמל יישוב                0
נפה                      0
מעמד מונציפאלי           0
דת יישוב                 0
סך הכל אוכלוסייה 2018    0
ערבים                    0
שנת ייסוד                0
צורת יישוב שוטפת         0
השתייכות ארגונית         0
קואורדינטות              0
גובה                     0
מרחב משטרה               0
dtype: int64

Data Types

In [44]:
cities.dtypes
Out[44]:
שם יישוב                  object
סמל יישוב                  int64
נפה                        int64
מעמד מונציפאלי            object
דת יישוב                  object
סך הכל אוכלוסייה 2018    float64
ערבים                    float64
שנת ייסוד                float64
צורת יישוב שוטפת           int64
השתייכות ארגונית         float64
קואורדינטות              float64
גובה                     float64
מרחב משטרה                object
dtype: object

The data types are changed to match the most suitable type.

In [45]:
cities = cities.astype('category')
cities = cities.astype({'סמל יישוב': 'str', 'סך הכל אוכלוסייה 2018': 'float64', 'ערבים': 'float64', 
                        'שנת ייסוד': 'int64', 'קואורדינטות': 'float64', 'גובה': 'float64'})
In [46]:
cities.dtypes
Out[46]:
שם יישוב                 category
סמל יישוב                  object
נפה                      category
מעמד מונציפאלי           category
דת יישוב                 category
סך הכל אוכלוסייה 2018     float64
ערבים                     float64
שנת ייסוד                   int64
צורת יישוב שוטפת         category
השתייכות ארגונית         category
קואורדינטות               float64
גובה                      float64
מרחב משטרה               category
dtype: object

1.3 Merging the Datasets

The crimes dataset does not include data for cities that belong to regional councils. Therefore, we will combine the data of all the cities that belong to the same regional council.

In [47]:
regional = cities.groupby(['מעמד מונציפאלי'])

For categorical values we decided to take the mode (i.e. most frequent value) for each regional council.

In [48]:
temp1 = regional[['נפה', 'דת יישוב', 'צורת יישוב שוטפת', 'השתייכות ארגונית', 'מרחב משטרה']].agg(lambda x:x.value_counts().index[0])

For variables that represent the number of citizens, we take the sum.

In [49]:
temp2 = regional[['סך הכל אוכלוסייה 2018', 'ערבים']].sum()

For any other numerical values, we take the mean.

In [50]:
temp3 = regional[['שנת ייסוד', 'קואורדינטות', 'גובה']].mean()

We then combine this data to a regional dataframe.

In [51]:
regional = pd.concat([temp1, temp2, temp3], axis=1)
regional = regional[1:55]
regional.head()
Out[51]:
נפה דת יישוב צורת יישוב שוטפת השתייכות ארגונית מרחב משטרה סך הכל אוכלוסייה 2018 ערבים שנת ייסוד קואורדינטות גובה
מעמד מונציפאלי
1.0 21 1 330 15.0 1.00046e+07 29290.600000 2562.60000 1944.257143 2.536552e+09 237.400000
2.0 21 1 310 19.0 1.00045e+07 23784.400000 8245.40000 1962.933333 2.421200e+09 486.406667
3.0 22 1 310 19.0 1.00046e+07 34320.093268 4778.22488 1950.440000 2.395178e+09 154.407602
4.0 24 1 330 15.0 1.00043e+07 88644.400000 30497.00000 1950.975000 2.138302e+09 92.760000
6.0 22 1 330 15.0 1.00046e+07 26500.693268 2910.82488 1939.814815 2.520133e+09 -99.200369

To get the names and corresponding IDs of each regional council, we use the mapping file and match the numbers to the names from the crimes dataset.

In [52]:
# get regional council names from the crimes dataset
names = crimes.loc[crimes['סמל יישוב'] == -2]['שם יישוב'].unique()
In [53]:
# get mapping of 'מעמד מוניציפלי'
mapping = pd.read_excel('mapping.xlsx', sheet_name='מעמד מוניציפלי')
In [54]:
# get the numbers of the regional councils
mapping = mapping.loc[11:66].sort_values(by=' מעמד מוניציפלי')
nums = mapping['Unnamed: 1'].values.astype('float64')
In [55]:
# create df
df = pd.DataFrame([names, nums, nums], index=['שם יישוב', 'מעמד מונציפאלי', 'סמל יישוב']).T
df.head()
Out[55]:
שם יישוב מעמד מונציפאלי סמל יישוב
0 מועצה אזורית אל-בטוף 69 69
1 מועצה אזורית אל קסום 65 65
2 מועצה אזורית אלונה 45 45
3 מועצה אזורית אשכול 38 38
4 מועצה אזורית באר טוביה 33 33
In [56]:
# merge names and numbers columns to regional dataset
regionals = pd.merge(regional, df, on='מעמד מונציפאלי')
regionals.head()
Out[56]:
מעמד מונציפאלי נפה דת יישוב צורת יישוב שוטפת השתייכות ארגונית מרחב משטרה סך הכל אוכלוסייה 2018 ערבים שנת ייסוד קואורדינטות גובה שם יישוב סמל יישוב
0 1 21 1 330 15.0 1.00046e+07 29290.600000 2562.60000 1944.257143 2.536552e+09 237.400000 מועצה אזורית הגליל העליון 1
1 2 21 1 310 19.0 1.00045e+07 23784.400000 8245.40000 1962.933333 2.421200e+09 486.406667 מועצה אזורית מרום הגליל 2
2 3 22 1 310 19.0 1.00046e+07 34320.093268 4778.22488 1950.440000 2.395178e+09 154.407602 מועצה אזורית הגליל התחתון 3
3 4 24 1 330 15.0 1.00043e+07 88644.400000 30497.00000 1950.975000 2.138302e+09 92.760000 מועצה אזורית מטה אשר 4
4 6 22 1 330 15.0 1.00046e+07 26500.693268 2910.82488 1939.814815 2.520133e+09 -99.200369 מועצה אזורית עמק הירדן 6

Next, we drop any cities that belong to reginal councils from the cities dataset and add the new rows with the regional data.

In [57]:
cities = cities.loc[(cities['מעמד מונציפאלי'] == 99) | (cities['מעמד מונציפאלי'] == 0) | (cities['מעמד מונציפאלי'] == 'חסר מעמד')]
In [58]:
cities.append(regionals, ignore_index=True)
Out[58]:
שם יישוב סמל יישוב נפה מעמד מונציפאלי דת יישוב סך הכל אוכלוסייה 2018 ערבים שנת ייסוד צורת יישוב שוטפת השתייכות ארגונית קואורדינטות גובה מרחב משטרה
0 אבו ג'ווייעד (שבט) 967 62 חסר מעמד 3 25830.400000 0.00000 1969.000000 460 19.0 2.040057e+09 86.000000 1.50037e+07
1 אבו גוש 472 11 99 2 7543.000000 7446.00000 1949.000000 280 19.0 2.105263e+09 598.000000 1.00025e+07
2 אבו סנאן 473 24 99 2 13915.000000 13887.00000 1968.000000 270 19.0 2.160776e+09 19.000000 1.00043e+07
3 אבו סריחאן (שבט) 935 62 חסר מעמד 3 16046.000000 0.00000 1943.000000 460 19.0 1.865057e+09 45.000000 1.00019e+07
4 אבו עבדון (שבט) 958 62 חסר מעמד 3 15094.200000 0.00000 1937.000000 460 19.0 1.850058e+09 70.000000 1.00019e+07
... ... ... ... ... ... ... ... ... ... ... ... ... ...
329 מועצה אזורית מטה בנימין 73 74 73 1 74418.893268 1550.02488 1982.344828 370 19.0 2.180452e+09 461.027243 1.00041e+07
330 מועצה אזורית מגילות ים המלח 74 75 74 1 8978.893268 1459.02488 1973.428571 330 15.0 2.373208e+09 -254.601422 1.00033e+07
331 מועצה אזורית עמק לוד 75 75 75 1 16485.293268 1515.42488 1974.000000 310 1.0 2.429336e+09 -112.600452 1.00033e+07
332 מועצה אזורית גוש עציון 76 76 76 1 31157.893268 1524.02488 1976.533333 370 19.0 2.155263e+09 719.652670 1.50006e+07
333 מועצה אזורית הר חברון 78 77 78 1 16988.893268 1460.82488 1981.000000 370 11.0 2.044042e+09 616.549378 1.50035e+07

334 rows × 13 columns

In [59]:
cities.describe()
Out[59]:
סך הכל אוכלוסייה 2018 ערבים שנת ייסוד קואורדינטות גובה
count 280.000000 280.000000 280.000000 2.800000e+02 280.000000
mean 31127.246377 6453.127939 1951.728571 2.078102e+09 153.832253
std 73018.767169 22942.032079 23.055595 2.211821e+08 191.184953
min 248.000000 0.000000 1870.000000 1.596762e+09 -227.000000
25% 5953.750000 10.000000 1944.500000 1.919671e+09 25.800000
50% 11987.000000 296.500000 1954.000000 2.072078e+09 89.500000
75% 25878.550000 6691.200000 1968.000000 2.227550e+09 216.700000
max 919438.000000 349572.000000 1998.000000 2.730079e+09 990.000000

We see that we have some outliers in the column סך הכל אוכלוסייה 2018 (the max value is significantly larger than the 75% value). Thus, we drop the rows that have more than 300000 citizens.

In [60]:
cities = cities.loc[cities['סך הכל אוכלוסייה 2018'] < 300000]

In the crimes dataset, we update the IDs of the regional councils.

In [61]:
# get regional councils crimes data
crimes_regional = crimes_train_wide.loc[crimes_train_wide['סמל יישוב'] == -2].sort_values(by='שם יישוב')
In [62]:
# update the ID
crimes_regional['מספר יישוב'] = nums
crimes_regional.drop(['סמל יישוב'], axis=1, inplace=True)
crimes_regional = crimes_regional.rename(columns={'מספר יישוב': 'סמל יישוב'})
In [63]:
crimes_regional.head()
Out[63]:
עבירה שם יישוב סעיפי הגדרה עבירות בטחון עבירות כלכליות עבירות כלפי המוסר עבירות כלפי הרכוש עבירות מין עבירות מנהליות עבירות מרמה עבירות נגד אדם עבירות נגד גוף עבירות סדר ציבורי עבירות רשוי עבירות תנועה שאר עבירות סמל יישוב
118 מועצה אזורית אל קסום 0.2 25.0 0.2 14.2 76.6 2.6 0.0 4.4 2.2 78.6 223.8 1.2 4.2 1.0 69.0
119 מועצה אזורית אל-בטוף 0.2 14.0 0.2 11.2 61.4 2.4 0.0 3.4 0.8 57.8 89.4 1.0 2.2 0.0 65.0
120 מועצה אזורית אלונה 0.2 0.0 0.0 4.4 61.6 0.8 0.0 0.8 0.0 4.8 17.0 0.0 0.0 0.0 45.0
121 מועצה אזורית אשכול 0.4 59.0 0.8 88.6 204.8 4.8 0.2 25.0 2.4 54.2 147.4 2.6 5.8 0.4 38.0
122 מועצה אזורית באר טוביה 2.4 8.2 0.6 79.6 705.0 12.8 0.0 32.4 1.6 123.0 300.6 7.2 3.8 0.4 33.0

Now we drop the previous regional councils data and append the new data.

In [64]:
# drop regional councils data
crimes_train_wide = crimes_train_wide.loc[crimes_train_wide['סמל יישוב'] != -2]
In [65]:
# append new data
crimes_train_wide.append(crimes_regional, ignore_index=True)
Out[65]:
שם יישוב סמל יישוב סעיפי הגדרה עבירות בטחון עבירות כלכליות עבירות כלפי המוסר עבירות כלפי הרכוש עבירות מין עבירות מנהליות עבירות מרמה עבירות נגד אדם עבירות נגד גוף עבירות סדר ציבורי עבירות רשוי עבירות תנועה שאר עבירות
0 אבו גוש 472.0 0.0 9.6 0.4 9.2 96.4 5.6 0.0 11.2 1.6 120.2 193.8 1.6 2.6 0.6
1 אבו סנאן 473.0 0.0 12.0 0.2 38.0 77.0 1.2 0.0 5.8 0.2 71.0 133.4 2.0 2.0 0.2
2 אבן יהודה 182.0 1.0 0.4 0.2 14.8 260.0 7.4 0.4 12.4 0.4 49.2 104.6 0.4 1.2 0.2
3 אום אל פחם 2710.0 1.2 118.4 4.2 112.0 433.4 7.0 0.4 47.0 7.0 240.6 655.6 5.6 11.6 0.6
4 אופקים 31.0 1.6 9.0 1.6 119.0 349.0 26.0 0.2 33.4 2.2 223.6 439.2 11.4 7.2 1.4
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
251 מועצה אזורית שדות נגב 39.0 0.4 2.6 0.4 40.6 181.8 4.4 0.0 8.0 0.2 42.8 106.4 1.2 3.8 0.8
252 מועצה אזורית שומרון 72.0 0.4 177.8 0.8 44.4 261.0 19.6 0.2 31.6 8.4 132.4 275.2 3.4 6.0 2.0
253 מועצה אזורית שער הנגב 37.0 0.4 34.6 1.8 25.0 130.4 2.2 0.0 10.0 1.8 32.2 81.8 4.2 1.6 0.6
254 מועצה אזורית שפיר 34.0 0.6 3.8 0.0 18.8 232.4 4.0 0.2 7.8 1.2 54.0 126.8 2.6 2.4 0.4
255 מועצה אזורית תמר 51.0 0.0 3.4 0.4 17.6 185.0 11.6 0.2 13.0 1.4 63.0 105.4 2.6 1.4 1.4

256 rows × 16 columns

One-Hot Encoding

We create dummy variables for all categorical columns.

In [66]:
# create dummy variables
cols = {'נפה', 'מעמד מונציפאלי', 'דת יישוב', 'צורת יישוב שוטפת', 'השתייכות ארגונית', 'מרחב משטרה'}
dummy_data = pd.get_dummies(cities, columns=cols, drop_first=True)
In [67]:
dummy_data.head()
Out[67]:
שם יישוב סמל יישוב סך הכל אוכלוסייה 2018 ערבים שנת ייסוד קואורדינטות גובה דת יישוב_2.0 דת יישוב_3.0 דת יישוב_4.0 ... נפה_53 נפה_61 נפה_62 נפה_71 נפה_72 נפה_73 נפה_74 נפה_75 נפה_76 נפה_77
0 אבו ג'ווייעד (שבט) 967 25830.4 0.0 1969 2.040057e+09 86.0 0 1 0 ... 0 0 1 0 0 0 0 0 0 0
1 אבו גוש 472 7543.0 7446.0 1949 2.105263e+09 598.0 1 0 0 ... 0 0 0 0 0 0 0 0 0 0
2 אבו סנאן 473 13915.0 13887.0 1968 2.160776e+09 19.0 1 0 0 ... 0 0 0 0 0 0 0 0 0 0
3 אבו סריחאן (שבט) 935 16046.0 0.0 1943 1.865057e+09 45.0 0 1 0 ... 0 0 1 0 0 0 0 0 0 0
4 אבו עבדון (שבט) 958 15094.2 0.0 1937 1.850058e+09 70.0 0 1 0 ... 0 0 1 0 0 0 0 0 0 0

5 rows × 203 columns

Merging

In order to merge the datasets, we make sure both IDs are the same type.

In [68]:
# change ID type to float
dummy_data['סמל יישוב'] = dummy_data['סמל יישוב'].astype('float64')

In addition, we drop the city name variable because it is not exactly the same in both datasets.

In [69]:
# drop city name
dummy_data.drop('שם יישוב', axis=1, inplace=True)
crimes_train_wide.drop('שם יישוב', axis=1, inplace=True)
crimes_test_wide.drop('שם יישוב', axis=1, inplace=True)

Lastly, we merge both datasets using the ID סמל יישוב.

In [70]:
train_df = pd.merge(dummy_data, crimes_train_wide, on='סמל יישוב')
In [71]:
train_df.head()
Out[71]:
סמל יישוב סך הכל אוכלוסייה 2018 ערבים שנת ייסוד קואורדינטות גובה דת יישוב_2.0 דת יישוב_3.0 דת יישוב_4.0 דת יישוב_אחר ... עבירות כלפי הרכוש עבירות מין עבירות מנהליות עבירות מרמה עבירות נגד אדם עבירות נגד גוף עבירות סדר ציבורי עבירות רשוי עבירות תנועה שאר עבירות
0 472.0 7543.0 7446.0 1949 2.105263e+09 598.0 1 0 0 0 ... 96.4 5.6 0.0 11.2 1.6 120.2 193.8 1.6 2.6 0.6
1 473.0 13915.0 13887.0 1968 2.160776e+09 19.0 1 0 0 0 ... 77.0 1.2 0.0 5.8 0.2 71.0 133.4 2.0 2.0 0.2
2 182.0 13700.0 2.0 1932 1.894469e+09 10.0 0 0 0 0 ... 260.0 7.4 0.4 12.4 0.4 49.2 104.6 0.4 1.2 0.2
3 2710.0 55182.0 55138.0 1950 2.145171e+09 139.0 1 0 0 0 ... 433.4 7.0 0.4 47.0 7.0 240.6 655.6 5.6 11.6 0.6
4 31.0 29021.0 201.0 1955 1.638258e+09 113.0 0 0 0 0 ... 349.0 26.0 0.2 33.4 2.2 223.6 439.2 11.4 7.2 1.4

5 rows × 216 columns

In [72]:
train_df.describe()
Out[72]:
סמל יישוב סך הכל אוכלוסייה 2018 ערבים שנת ייסוד קואורדינטות גובה דת יישוב_2.0 דת יישוב_3.0 דת יישוב_4.0 דת יישוב_אחר ... עבירות כלפי הרכוש עבירות מין עבירות מנהליות עבירות מרמה עבירות נגד אדם עבירות נגד גוף עבירות סדר ציבורי עבירות רשוי עבירות תנועה שאר עבירות
count 198.000000 198.000000 198.000000 198.000000 1.980000e+02 198.000000 198.000000 198.0 198.000000 198.0 ... 198.000000 198.000000 198.000000 198.000000 198.000000 198.000000 198.000000 198.000000 198.000000 198.000000
mean 2618.595960 33341.762626 6882.914141 1951.116162 2.099243e+09 151.000000 0.409091 0.0 0.030303 0.0 ... 491.743434 20.941414 0.390909 58.995960 2.244444 221.424242 440.804040 4.966667 5.141414 0.511111
std 2857.555147 48768.062404 11535.590354 25.124495 2.359100e+08 207.403552 0.492912 0.0 0.171854 0.0 ... 968.807124 37.990400 1.397648 123.879486 3.279105 379.704146 721.604492 8.621982 7.948221 0.969681
min 26.000000 1241.000000 0.000000 1878.000000 1.596762e+09 -227.000000 0.000000 0.0 0.000000 0.0 ... 3.200000 0.000000 0.000000 0.200000 0.000000 2.200000 6.800000 0.000000 0.000000 0.000000
25% 517.250000 7602.000000 45.000000 1942.000000 1.920846e+09 20.000000 0.000000 0.0 0.000000 0.0 ... 60.150000 2.200000 0.000000 5.800000 0.400000 35.450000 77.100000 0.800000 1.000000 0.000000
50% 1176.500000 17215.000000 524.500000 1953.000000 2.088021e+09 62.000000 0.000000 0.0 0.000000 0.0 ... 139.900000 5.300000 0.000000 15.000000 1.000000 84.400000 181.500000 2.200000 2.400000 0.200000
75% 3745.000000 36576.750000 10330.250000 1969.000000 2.260899e+09 228.750000 1.000000 0.0 0.000000 0.0 ... 420.700000 22.650000 0.200000 51.500000 2.600000 239.100000 466.350000 5.000000 5.850000 0.600000
max 9800.000000 283640.000000 76918.000000 1998.000000 2.730079e+09 990.000000 1.000000 0.0 1.000000 0.0 ... 6662.000000 266.200000 12.000000 1093.400000 25.600000 2490.200000 4986.600000 68.600000 71.800000 8.400000

8 rows × 216 columns

In [73]:
# merge also test dataset
test_df = pd.merge(dummy_data, crimes_test_wide, on='סמל יישוב')

We make another merged dataset that does not include dummy variables which will be used in some specific cases later.

In [74]:
cities['סמל יישוב'] = cities['סמל יישוב'].astype('float64')
cities.drop('שם יישוב', axis=1, inplace=True)
df_clustering = pd.merge(cities, crimes_train_wide, on='סמל יישוב')

Scaling

We standardize the columns for better performance of the algorithms.

In [75]:
train_features = train_df.loc[:, 'סך הכל אוכלוסייה 2018': 'שאר עבירות']
ID = train_df.loc[:, 'סמל יישוב']
In [76]:
# scaling train dataset
scaler = StandardScaler()
scaled = pd.DataFrame(scaler.fit_transform(train_features), columns = train_features.columns)
In [77]:
train_df_scaled = pd.concat([ID, scaled], axis=1)
In [78]:
train_df_scaled.head()
Out[78]:
סמל יישוב סך הכל אוכלוסייה 2018 ערבים שנת ייסוד קואורדינטות גובה דת יישוב_2.0 דת יישוב_3.0 דת יישוב_4.0 דת יישוב_אחר ... עבירות כלפי הרכוש עבירות מין עבירות מנהליות עבירות מרמה עבירות נגד אדם עבירות נגד גוף עבירות סדר ציבורי עבירות רשוי עבירות תנועה שאר עבירות
0 472.0 -0.530350 0.048937 -0.084441 0.025586 2.160682 1.20185 0.0 -0.176777 0.0 ... -0.409107 -0.404847 -0.280400 -0.386804 -0.197029 -0.267263 -0.343166 -0.391465 -0.320557 0.091901
1 473.0 -0.399360 0.608711 0.673711 0.261496 -0.638054 1.20185 0.0 -0.176777 0.0 ... -0.429182 -0.520959 -0.280400 -0.430506 -0.625057 -0.397166 -0.427081 -0.344954 -0.396237 -0.321652
2 182.0 -0.403780 -0.598006 -0.762786 -0.870218 -0.681557 -0.83205 0.0 -0.176777 0.0 ... -0.239811 -0.357347 0.006521 -0.377093 -0.563910 -0.454725 -0.467093 -0.530997 -0.497143 -0.321652
3 2710.0 0.448974 4.193752 -0.044538 0.195180 -0.058005 1.20185 0.0 -0.176777 0.0 ... -0.060375 -0.367902 0.006521 -0.097081 1.453936 0.050630 0.298419 0.073642 0.814642 0.091901
4 31.0 -0.088823 -0.580712 0.154976 -1.959024 -0.183682 -0.83205 0.0 -0.176777 0.0 ... -0.147713 0.133492 -0.136939 -0.207144 -0.013588 0.005745 -0.002229 0.748046 0.259656 0.919005

5 rows × 216 columns

In [79]:
# scale test df
test_features = test_df.loc[:, 'סך הכל אוכלוסייה 2018': 'שאר עבירות']
ID = test_df.loc[:, 'סמל יישוב']

# scaling train dataset
scaler = StandardScaler()
scaled = pd.DataFrame(scaler.fit_transform(test_features), columns = test_features.columns)

test_df_scaled = pd.concat([ID, scaled], axis=1)

2 Data Visualization

The following heatmap represents correlation of features. It seems like some types of crimes are very correlated.

In [80]:
temp = df_clustering.drop(['סמל יישוב'], axis=1)
temp_val=temp.corr().columns.values 
temp_val[0]='סך הכל אוכלוסייה 8102'# to make sure that the reverse will work also on the number 2018
plt.figure(figsize=(15, 12))
sns.heatmap(temp.corr(), cmap='Blues', xticklabels=[i[::-1] for i in temp_val],
            yticklabels=[i[::-1] for i in temp_val])
plt.title('Heatmap of Cities Dataset')
plt.show()

The following plot shows the distribution of the סך הכל אוכלוסייה 2018 feature.

In [81]:
plt.figure(figsize=(10, 8))
sns.histplot(df_clustering, x='סך הכל אוכלוסייה 2018', color='#b3cde3')
plt.grid(axis='y', alpha=0.5, color='lightgrey')
plt.title('Histogram of Population')
plt.xlabel('Population')
plt.show()

The following plot shows the overall number of crimes for the 10 cities with the most crimes.

In [82]:
top10_df = crimes.groupby('שם יישוב').sum()[[2019]].sort_values(by=[2019], ascending=False).head(11)
top10_df = top10_df.drop('אחר')
plt.figure(figsize=(10, 8))
plt.bar([i[::-1] for i in top10_df.index], top10_df[2019], color='#b3cde3', edgecolor='black')
plt.title('Cities with Most Overall Crimes in 2019')
plt.xlabel('City')
plt.ylabel('Overall Number of Crimes')
plt.grid(axis='y', alpha=0.5, color='lightgrey')
plt.show()

3 Clustering

In this section, we attempt to identify clusters in the data, using K-Means and GMM.

We define to functions to avoid code repetition:

The fit_model function recieves as input the data, the clustering model ('kmeans' or 'gmm') and the number of desired clusters.

The function fits the model to the data and returns the data with an added column of clusters.

In [83]:
def fit_model(X, mod, k):
    
    if mod == 'kmeans':
        model = KMeans(n_clusters=k, random_state=RSEED)
    elif mod == 'gmm':
        model = GaussianMixture(n_components=k, random_state=RSEED)
    else:
        raise ValueError('Model not found')
        
    model.fit(X)
    clusters = model.predict(X)
    X["Cluster"] = clusters
    
    return X

The perform_pca function, recieves as input the data and performs PCA with 2 components. It returns the dataframe with two columns with the first and second principal component. This is used to visualize the clusters.

In [84]:
def perform_pca(X):
    pca_2d = PCA(n_components=2, random_state=RSEED)
    PCs_2d = pd.DataFrame(pca_2d.fit_transform(X.drop(["Cluster"], axis=1)))
    PCs_2d.columns = ["PC1_2d", "PC2_2d"]
    X = pd.concat([X, PCs_2d], axis=1, join='inner')
    
    return X

3.1 K-Means

The first clustering algorithm that we implement is K-Means. We use $K=3,4,5$.

K = 3

In [85]:
df = fit_model(train_df.drop(['סמל יישוב'], axis=1), 'kmeans', 3)
X = perform_pca(df)
kmeans_df3 = X.copy()
In [86]:
cluster0 = X[X["Cluster"] == 0]
cluster1 = X[X["Cluster"] == 1]
cluster2 = X[X["Cluster"] == 2]
In [87]:
trace0 = go.Scatter(x = cluster0["PC1_2d"], y = cluster0["PC2_2d"],
                    mode = "markers", name = "Cluster 0",
                    marker = dict(color = '#fbb4ae'), text = None)

trace1 = go.Scatter(x = cluster1["PC1_2d"], y = cluster1["PC2_2d"],
                    mode = "markers", name = "Cluster 1",
                    marker = dict(color = '#b3cde3'), text = None)

trace2 = go.Scatter(x = cluster2["PC1_2d"], y = cluster2["PC2_2d"],
                    mode = "markers", name = "Cluster 2",
                    marker = dict(color = '#ccebc5'), text = None)

data = [trace0, trace1, trace2]

title = "Visualizing K-Means with K=3 Using PCA"

layout = dict(title = title,
              xaxis = dict(title='PC1', ticklen=5, zeroline=False),
              yaxis = dict(title='PC2', ticklen=5, zeroline=False),
              plot_bgcolor = '#f0f0f0')

fig = dict(data = data, layout = layout)

iplot(fig)

K = 4

In [88]:
df = fit_model(train_df.drop(['סמל יישוב'], axis=1), 'kmeans', 4)
X = perform_pca(df)
kmeans_df4 = X.copy()
In [89]:
cluster0 = X[X["Cluster"] == 0]
cluster1 = X[X["Cluster"] == 1]
cluster2 = X[X["Cluster"] == 2]
cluster3 = X[X["Cluster"] == 3]
In [90]:
trace0 = go.Scatter(x = cluster0["PC1_2d"], y = cluster0["PC2_2d"],
                    mode = "markers", name = "Cluster 0",
                    marker = dict(color = '#fbb4ae'), text = None)

trace1 = go.Scatter(x = cluster1["PC1_2d"], y = cluster1["PC2_2d"],
                    mode = "markers", name = "Cluster 1",
                    marker = dict(color = '#b3cde3'), text = None)

trace2 = go.Scatter(x = cluster2["PC1_2d"], y = cluster2["PC2_2d"],
                    mode = "markers", name = "Cluster 2",
                    marker = dict(color = '#ccebc5'), text = None)

trace3 = go.Scatter(x = cluster3["PC1_2d"], y = cluster3["PC2_2d"],
                    mode = "markers", name = "Cluster 3",
                    marker = dict(color = '#decbe4'), text = None)

data = [trace0, trace1, trace2, trace3]

title = "Visualizing K-Means with K=4 Using PCA"

layout = dict(title = title,
              xaxis = dict(title='PC1', ticklen=5, zeroline=False),
              yaxis = dict(title='PC2', ticklen=5, zeroline=False),
              plot_bgcolor = '#f0f0f0')

fig = dict(data = data, layout = layout)

iplot(fig)

K = 5

In [91]:
df = fit_model(train_df.drop(['סמל יישוב'], axis=1), 'kmeans', 5)
X = perform_pca(df)
kmeans_df5 = X.copy()
In [92]:
cluster0 = X[X["Cluster"] == 0]
cluster1 = X[X["Cluster"] == 1]
cluster2 = X[X["Cluster"] == 2]
cluster3 = X[X["Cluster"] == 3]
cluster4 = X[X["Cluster"] == 4]
In [93]:
trace0 = go.Scatter(x = cluster0["PC1_2d"], y = cluster0["PC2_2d"],
                    mode = "markers", name = "Cluster 0",
                    marker = dict(color = '#fbb4ae'), text = None)

trace1 = go.Scatter(x = cluster1["PC1_2d"], y = cluster1["PC2_2d"],
                    mode = "markers", name = "Cluster 1",
                    marker = dict(color = '#b3cde3'), text = None)

trace2 = go.Scatter(x = cluster2["PC1_2d"], y = cluster2["PC2_2d"],
                    mode = "markers", name = "Cluster 2",
                    marker = dict(color = '#ccebc5'), text = None)

trace3 = go.Scatter(x = cluster3["PC1_2d"], y = cluster3["PC2_2d"],
                    mode = "markers", name = "Cluster 3",
                    marker = dict(color = '#decbe4'), text = None)

trace4 = go.Scatter(x = cluster4["PC1_2d"], y = cluster4["PC2_2d"],
                    mode = "markers", name = "Cluster 4",
                    marker = dict(color = '#fed9a6'), text = None)

data = [trace0, trace1, trace2, trace3, trace4]

title = "Visualizing K-Means with K=5 Using PCA"

layout = dict(title = title,
              xaxis = dict(title='PC1', ticklen=5, zeroline=False),
              yaxis = dict(title='PC2', ticklen=5, zeroline=False),
              plot_bgcolor = '#f0f0f0')

fig = dict(data = data, layout = layout)

iplot(fig)

3.2 GMM

The second clustering algorithm that we implement is GMM. Again, we use $K=3,4,5$.

K = 3

In [94]:
df = fit_model(train_df.drop(['סמל יישוב'], axis=1), 'gmm', 3)
X = perform_pca(df)
gmm_df3 = X.copy()
In [95]:
cluster0 = X[X["Cluster"] == 0]
cluster1 = X[X["Cluster"] == 1]
cluster2 = X[X["Cluster"] == 2]
In [96]:
trace0 = go.Scatter(x = cluster0["PC1_2d"], y = cluster0["PC2_2d"],
                    mode = "markers", name = "Cluster 0",
                    marker = dict(color = '#fbb4ae'), text = None)

trace1 = go.Scatter(x = cluster1["PC1_2d"], y = cluster1["PC2_2d"],
                    mode = "markers", name = "Cluster 1",
                    marker = dict(color = '#b3cde3'), text = None)

trace2 = go.Scatter(x = cluster2["PC1_2d"], y = cluster2["PC2_2d"],
                    mode = "markers", name = "Cluster 2",
                    marker = dict(color = '#ccebc5'), text = None)

data = [trace0, trace1, trace2]

title = "Visualizing GMM with K=3 Using PCA"

layout = dict(title = title,
              xaxis = dict(title='PC1', ticklen=5, zeroline=False),
              yaxis = dict(title='PC2', ticklen=5, zeroline=False),
              plot_bgcolor = '#f0f0f0')

fig = dict(data = data, layout = layout)

iplot(fig)

K = 4

In [97]:
df = fit_model(train_df.drop(['סמל יישוב'], axis=1), 'gmm', 4)
X = perform_pca(df)
gmm_df4 = X.copy()
In [98]:
cluster0 = X[X["Cluster"] == 0]
cluster1 = X[X["Cluster"] == 1]
cluster2 = X[X["Cluster"] == 2]
cluster3 = X[X["Cluster"] == 3]
In [99]:
trace0 = go.Scatter(x = cluster0["PC1_2d"], y = cluster0["PC2_2d"],
                    mode = "markers", name = "Cluster 0",
                    marker = dict(color = '#fbb4ae'), text = None)

trace1 = go.Scatter(x = cluster1["PC1_2d"], y = cluster1["PC2_2d"],
                    mode = "markers", name = "Cluster 1",
                    marker = dict(color = '#b3cde3'), text = None)

trace2 = go.Scatter(x = cluster2["PC1_2d"], y = cluster2["PC2_2d"],
                    mode = "markers", name = "Cluster 2",
                    marker = dict(color = '#ccebc5'), text = None)

trace3 = go.Scatter(x = cluster3["PC1_2d"], y = cluster3["PC2_2d"],
                    mode = "markers", name = "Cluster 3",
                    marker = dict(color = '#decbe4'), text = None)

data = [trace0, trace1, trace2, trace3]

title = "Visualizing GMM with K=4 Using PCA"

layout = dict(title = title,
              xaxis = dict(title='PC1', ticklen=5, zeroline=False),
              yaxis = dict(title='PC2', ticklen=5, zeroline=False),
              plot_bgcolor = '#f0f0f0')

fig = dict(data = data, layout = layout)

iplot(fig)

K = 5

In [100]:
df = fit_model(train_df.drop(['סמל יישוב'], axis=1), 'gmm', 5)
X = perform_pca(df)
gmm_df5 = X.copy()
In [101]:
cluster0 = X[X["Cluster"] == 0]
cluster1 = X[X["Cluster"] == 1]
cluster2 = X[X["Cluster"] == 2]
cluster3 = X[X["Cluster"] == 3]
cluster4 = X[X["Cluster"] == 4]
In [102]:
trace0 = go.Scatter(x = cluster0["PC1_2d"], y = cluster0["PC2_2d"],
                    mode = "markers", name = "Cluster 0",
                    marker = dict(color = '#fbb4ae'), text = None)

trace1 = go.Scatter(x = cluster1["PC1_2d"], y = cluster1["PC2_2d"],
                    mode = "markers", name = "Cluster 1",
                    marker = dict(color = '#b3cde3'), text = None)

trace2 = go.Scatter(x = cluster2["PC1_2d"], y = cluster2["PC2_2d"],
                    mode = "markers", name = "Cluster 2",
                    marker = dict(color = '#ccebc5'), text = None)

trace3 = go.Scatter(x = cluster3["PC1_2d"], y = cluster3["PC2_2d"],
                    mode = "markers", name = "Cluster 3",
                    marker = dict(color = '#decbe4'), text = None)

trace4 = go.Scatter(x = cluster4["PC1_2d"], y = cluster4["PC2_2d"],
                    mode = "markers", name = "Cluster 4",
                    marker = dict(color = '#fed9a6'), text = None)

data = [trace0, trace1, trace2, trace3, trace4]

title = "Visualizing GMM with K=5 Using PCA"

layout = dict(title = title,
              xaxis = dict(title='PC1', ticklen=5, zeroline=False),
              yaxis = dict(title='PC2', ticklen=5, zeroline=False),
              plot_bgcolor = '#f0f0f0')

fig = dict(data = data, layout = layout)

iplot(fig)

3.3 Comparison and Conclusions

Both algorithms had very similar results (i.e. similar clusters).

For each cluster, we calculate the percentage of arabs from the overall population, and check whether it affects the percentage of עבירות נגד אדם and עבירות בטחון from the overall crimes.

We define the function create_res_df that creates a dataframe with the required percentages per cluster.

In [103]:
def create_res_df(df):
    dat = pd.concat([df_clustering, df['Cluster']], axis=1)
    # group df
    grouped = dat.groupby('Cluster')
    # calculate arab percentage
    percentage = grouped['ערבים'].sum() / grouped['סך הכל אוכלוסייה 2018'].sum()
    # calculate sum of all crimes
    cols = [col for col in df_clustering.columns if 'עבירות' in col]
    offenses = df_clustering[cols]
    crimes = grouped[offenses.columns].sum().sum(axis=1)
    # calculate crimes percentage
    crime1 = grouped['עבירות נגד אדם'].sum() / crimes
    crime2 = grouped['עבירות בטחון'].sum() / crimes
    
    # create res df
    res = pd.concat([percentage, crime1, crime2], axis=1).sort_values(by=0)
    res = res.rename(columns={0: 'אחוז הערבים מכלל האוכלוסייה', 1: 'אחוז עבירות נגד אדם מסך העבירות', 2: 'אחוז עבירות ביטחון מסך העבירות'})
    
    return res

K-Means

In [104]:
create_res_df(kmeans_df3)
Out[104]:
אחוז הערבים מכלל האוכלוסייה אחוז עבירות נגד אדם מסך העבירות אחוז עבירות ביטחון מסך העבירות
Cluster
1 0.073253 0.001421 0.008710
2 0.395321 0.001890 0.019777
0 0.544700 0.002276 0.025592
In [105]:
create_res_df(kmeans_df4)
Out[105]:
אחוז הערבים מכלל האוכלוסייה אחוז עבירות נגד אדם מסך העבירות אחוז עבירות ביטחון מסך העבירות
Cluster
2 0.065429 0.001368 0.008072
1 0.193534 0.001765 0.014376
0 0.311212 0.001899 0.017899
3 0.677749 0.002399 0.031196
In [106]:
create_res_df(kmeans_df5)
Out[106]:
אחוז הערבים מכלל האוכלוסייה אחוז עבירות נגד אדם מסך העבירות אחוז עבירות ביטחון מסך העבירות
Cluster
0 0.049249 0.001386 0.006648
2 0.174261 0.001576 0.013216
3 0.193534 0.001765 0.014376
4 0.410290 0.001971 0.021669
1 0.686127 0.002378 0.032119

GMM

In [107]:
create_res_df(gmm_df3)
Out[107]:
אחוז הערבים מכלל האוכלוסייה אחוז עבירות נגד אדם מסך העבירות אחוז עבירות ביטחון מסך העבירות
Cluster
0 0.073406 0.001423 0.008707
2 0.376249 0.001882 0.018948
1 0.560401 0.002212 0.026401
In [108]:
create_res_df(gmm_df4)
Out[108]:
אחוז הערבים מכלל האוכלוסייה אחוז עבירות נגד אדם מסך העבירות אחוז עבירות ביטחון מסך העבירות
Cluster
0 0.055309 0.001343 0.007296
3 0.193534 0.001765 0.014376
2 0.308182 0.001911 0.018780
1 0.677749 0.002399 0.031196
In [109]:
create_res_df(gmm_df5)
Out[109]:
אחוז הערבים מכלל האוכלוסייה אחוז עבירות נגד אדם מסך העבירות אחוז עבירות ביטחון מסך העבירות
Cluster
0 0.049249 0.001386 0.006648
2 0.174261 0.001576 0.013216
3 0.193534 0.001765 0.014376
4 0.410290 0.001971 0.021669
1 0.686127 0.002378 0.032119

We see that both algorithms had almost identical results (i.e. clusters) and in all of them we see that the higher the percentage of arabs from the overall population, the more crimes of type נגד אדם and ביטחון we have.

In our opinion, the best result is obtained by the K-Means algorithm with $K=3$ because it gives the best and clearest separation. It is not necessary to use more clusters because they do not lead to any improvment.


4 Random Forest

In this section, we attempt to predict the second most frequent crime for certain cities, using the Random Forest algorithm.

First we create a dataframe that includes only columns with the word עבירות.

In [110]:
cols = [col for col in train_df_scaled.columns if 'עבירות' in col]
offenses = train_df_scaled[cols]
offenses.head()
Out[110]:
עבירות בטחון עבירות כלכליות עבירות כלפי המוסר עבירות כלפי הרכוש עבירות מין עבירות מנהליות עבירות מרמה עבירות נגד אדם עבירות נגד גוף עבירות סדר ציבורי עבירות רשוי עבירות תנועה שאר עבירות
0 -0.327640 -0.403828 -0.566207 -0.409107 -0.404847 -0.280400 -0.386804 -0.197029 -0.267263 -0.343166 -0.391465 -0.320557 0.091901
1 -0.236044 -0.437452 -0.414116 -0.429182 -0.520959 -0.280400 -0.430506 -0.625057 -0.397166 -0.427081 -0.344954 -0.396237 -0.321652
2 -0.678757 -0.437452 -0.536634 -0.239811 -0.357347 0.006521 -0.377093 -0.563910 -0.454725 -0.467093 -0.530997 -0.497143 -0.321652
3 3.824702 0.235029 -0.023327 -0.060375 -0.367902 0.006521 -0.097081 1.453936 0.050630 0.298419 0.073642 0.814642 0.091901
4 -0.350539 -0.202084 0.013640 -0.147713 0.133492 -0.136939 -0.207144 -0.013588 0.005745 -0.002229 0.748046 0.259656 0.919005

Next, we select the required rows from the test dataset.

In [111]:
cities_sign = ([1139, 7400, 507, 2640, 43, 3780])
rf_test = test_df[test_df['סמל יישוב'].isin(cities_sign)]
rf_test
Out[111]:
סמל יישוב סך הכל אוכלוסייה 2018 ערבים שנת ייסוד קואורדינטות גובה דת יישוב_2.0 דת יישוב_3.0 דת יישוב_4.0 דת יישוב_אחר ... עבירות כלפי הרכוש עבירות מין עבירות מנהליות עבירות מרמה עבירות נגד אדם עבירות נגד גוף עבירות סדר ציבורי עבירות רשוי עבירות תנועה שאר עבירות
32 3780.0 56746.0 2.0 1985 2.107362e+09 568.0 0 0 0 0 ... 248.0 48.0 0.0 28.0 1.0 133.0 230.0 3.0 0.0 1.0
94 507.0 10029.0 10008.0 1968 2.155976e+09 34.0 1 0 0 0 ... 111.0 8.0 0.0 16.0 3.0 72.0 183.0 2.0 1.0 0.0
104 1139.0 46124.0 1507.0 1964 2.276276e+09 69.0 0 0 0 0 ... 474.0 47.0 0.0 96.0 3.0 285.0 563.0 15.0 9.0 1.0
118 43.0 1599.0 65.0 1896 2.541580e+09 260.0 0 0 0 0 ... 7.0 1.0 0.0 2.0 0.0 5.0 9.0 0.0 1.0 0.0
135 7400.0 217244.0 535.0 1929 1.872569e+09 -2.0 0 0 0 0 ... 3379.0 108.0 1.0 380.0 6.0 1366.0 2354.0 18.0 26.0 0.0
177 2640.0 56344.0 31.0 1950 1.961967e+09 29.0 0 0 0 0 ... 614.0 56.0 1.0 77.0 2.0 294.0 572.0 6.0 4.0 0.0

6 rows × 216 columns

Random Forest Regressor

We loop through the different crime types and using the train data, we predict the number of the specific crime type in 2019.

In [112]:
parameters_grid = {
    'max_depth': [50, 100, 150, 200],
    'min_samples_leaf': [1, 2],
    'min_samples_split': [2, 3],
    'n_estimators': [50, 100, 150],
    'max_features': [5, 10, 15, 20]
}

pred_res = []
MSE = []

# create a loop for each 'עבירות' that is in the dataframe
for i in range(len(offenses.columns)):
    
    # exclude the ith column from the train df
    train_cols = [col for col in train_df_scaled.columns if col not in [offenses.columns[i]]]
    # train df without the ith column
    temp_train = train_df_scaled[train_cols].drop(['סמל יישוב'], axis=1)
    
    # exclude the ith column from the test df
    test_cols = [col for col in test_df_scaled.columns if col not in [offenses.columns[i]]]
    # test df without the ith column
    temp_test = test_df_scaled[test_cols].drop(['סמל יישוב'], axis=1)
    
    # define model
    rf = RandomForestRegressor(random_state=RSEED)
    # define grid search
    grid_search = GridSearchCV(estimator=rf, param_grid=parameters_grid, cv=10, n_jobs=-1)
    # fit the model
    grid_search.fit(temp_train, train_df[offenses.columns[i]])
    # get best model
    best = grid_search.best_estimator_
    # predict using best model
    y_pred = best.predict(temp_test)
    
    # add predictions for the required cities
    pred_res.append(y_pred[rf_test.index])
    
    # add MSE of 
    MSE.append(mean_squared_error(test_df[offenses.columns[i]], y_pred))

The following table presents the predicted crime numbers for each test city.

In [113]:
rf_res = pd.DataFrame()
rf_res = rf_res.append(pred_res).T
rf_res.columns = offenses.columns
rf_res.index = ['כרמיאל','נתניה','כפר יאסיף','ראש העין','מטולה','ביתר עילית']
rf_res
Out[113]:
עבירות בטחון עבירות כלכליות עבירות כלפי המוסר עבירות כלפי הרכוש עבירות מין עבירות מנהליות עבירות מרמה עבירות נגד אדם עבירות נגד גוף עבירות סדר ציבורי עבירות רשוי עבירות תנועה שאר עבירות
כרמיאל 23.758 2.461809 112.504 327.408 27.882 0.405123 36.068 2.712 192.716 364.373190 4.504092 4.578667 0.471444
נתניה 19.848 1.168379 54.292 148.632 4.430 0.184360 18.252 1.152 105.168 207.544333 2.462743 3.178667 0.149889
כפר יאסיף 21.206 5.685175 238.418 613.972 37.552 0.356284 81.224 3.664 346.656 706.490622 10.868711 7.494667 0.573083
ראש העין 4.286 0.266638 13.828 48.624 1.400 0.161077 3.384 0.276 14.256 31.663489 0.602289 0.932000 0.204098
מטולה 35.088 14.943547 521.170 3094.276 117.950 1.393799 357.442 8.436 1046.744 2203.454552 22.399156 20.896000 2.450044
ביתר עילית 12.542 4.900271 167.296 780.524 32.820 0.412648 87.420 2.484 283.748 622.574993 6.876000 8.237333 0.558578

Now, we select the second most frequent crime according to the predictions and compare it to the actual data of 2019.

In [114]:
pred_second = rf_res.apply(lambda row: row.nlargest(2).values[-1], axis=1)
pred_ind = rf_res.apply(lambda row: row.nlargest(2).idxmin(), axis=1)
In [115]:
actual_second = rf_test[cols].apply(lambda row: row.nlargest(2).values[-1], axis=1)
actual_ind = rf_test[cols].apply(lambda row: row.nlargest(2).idxmin(), axis=1)
In [116]:
res_df = pd.DataFrame(pred_ind, columns=['עבירה חזויה'])
res_df['עבירה אמיתית'] = actual_ind.values
res_df['מספר העבירות החזוי'] = pred_second
res_df['מספר העבירות האמיתי'] = actual_second.values
In [117]:
res_df
Out[117]:
עבירה חזויה עבירה אמיתית מספר העבירות החזוי מספר העבירות האמיתי
כרמיאל עבירות כלפי הרכוש עבירות סדר ציבורי 327.408000 230.0
נתניה עבירות כלפי הרכוש עבירות כלפי המוסר 148.632000 121.0
כפר יאסיף עבירות כלפי הרכוש עבירות כלפי הרכוש 613.972000 474.0
ראש העין עבירות סדר ציבורי עבירות סדר ציבורי 31.663489 9.0
מטולה עבירות סדר ציבורי עבירות סדר ציבורי 2203.454552 2354.0
ביתר עילית עבירות סדר ציבורי עבירות סדר ציבורי 622.574993 572.0

We see that 5 out of the 6 predictions were accurate, which is a good result.

In addition, we calculated the mean MSE presented below.

In [118]:
rf_mse = pd.DataFrame(round(np.mean(MSE), 3), columns=['Mean MSE'], index=['Random Forest Regressor'])
rf_mse
Out[118]:
Mean MSE
Random Forest Regressor 4711.342

5 Ada-Boost

In this section, we attempt to predict the overall number of crimes in 2019 for certain cities.

First we create a dataframe that includes only columns with the word עבירות (both for the scaled and non-scaled data).

In [119]:
cols = [col for col in train_df_scaled.columns if 'עבירות' in col]
offenses = train_df_scaled[cols]
In [120]:
cols = [col for col in train_df.columns if 'עבירות' in col]
offenses_og = train_df[cols]

Next, we select the required rows from the test dataset.

In [121]:
cities_sign = ([2600, 9000, 7500, 6800, 26])
ab_test = test_df[test_df['סמל יישוב'].isin(cities_sign)]
ab_test
Out[121]:
סמל יישוב סך הכל אוכלוסייה 2018 ערבים שנת ייסוד קואורדינטות גובה דת יישוב_2.0 דת יישוב_3.0 דת יישוב_4.0 דת יישוב_אחר ... עבירות כלפי הרכוש עבירות מין עבירות מנהליות עבירות מרמה עבירות נגד אדם עבירות נגד גוף עבירות סדר ציבורי עבירות רשוי עבירות תנועה שאר עבירות
9 2600.0 51935.0 2260.0 1951 1.941739e+09 -5.0 0 0 0 0 ... 1816.0 122.0 5.0 179.0 8.0 1239.0 1824.0 21.0 27.0 2.0
22 9000.0 209002.0 5380.0 1948 1.796357e+09 196.0 0 0 0 0 ... 6499.0 364.0 8.0 535.0 25.0 2320.0 4612.0 57.0 73.0 1.0
137 7500.0 31057.0 31031.0 1958 2.283275e+09 178.0 1 0 0 0 ... 241.0 3.0 0.0 31.0 0.0 201.0 366.0 1.0 5.0 0.0
165 6800.0 58267.0 113.0 1925 2.104675e+09 1.0 0 0 0 0 ... 834.0 42.0 0.0 104.0 1.0 502.0 810.0 5.0 4.0 0.0
178 26.0 3120.0 39.0 1882 2.510676e+09 150.0 0 0 0 0 ... 72.0 2.0 0.0 19.0 2.0 22.0 58.0 1.0 8.0 0.0

5 rows × 216 columns

In addition, we define train and test datasets without the crimes columns.

In [122]:
cols = [col for col in train_df_scaled.columns if 'עבירות' not in col]
ab_train_scaled = train_df_scaled[cols]
In [123]:
cols = [col for col in test_df_scaled.columns if 'עבירות' not in col]
ab_test_scaled = test_df_scaled[cols]

The features of the train data do not include the test cities and thus they are removed, and the label is the sum of the crimes.

In [124]:
ab_train_X = ab_train_scaled.loc[~ab_train_scaled.index.isin(ab_test.index)]
ab_train_y = offenses_og.loc[~offenses_og.index.isin(ab_test.index)].sum(axis=1)

Ada-Boost Regressor

In [125]:
parameters_grid = {
    'n_estimators': [10, 20, 50, 100, 150, 200],
    'learning_rate': [0.0001, 0.001, 0.01, 0.1, 1.0, 1.5, 2]
}

ab = AdaBoostRegressor(random_state=RSEED)
grid_search = GridSearchCV(estimator=ab, param_grid=parameters_grid, cv=10, n_jobs=-1)

grid_search.fit(ab_train_X, ab_train_y)
best = grid_search.best_estimator_

y_pred = best.predict(ab_test_scaled.loc[ab_test_scaled.index.isin(ab_test.index)])

We now compare the predicted vs. the true overall sum of the crimes for each test city.

In [126]:
cols = [col for col in train_df.columns if 'עבירות' in col]
offenses_test = test_df[cols]
y_actual = offenses_test.loc[offenses_test.index.isin(ab_test.index)].sum(axis=1)
In [127]:
ab_res = pd.DataFrame(y_pred, columns={'מספר העבירות החזוי'})
ab_res.index=['אילת','באר שבע','סח\'נין','קרית אתא','ראש פינה']
ab_res['מספר העבירות האמיתי'] = offenses_test.sum(axis=1)[ab_test.index].values
ab_res
Out[127]:
מספר העבירות החזוי מספר העבירות האמיתי
אילת 2168.470588 5942.0
באר שבע 9685.250000 15286.0
סח'נין 1373.906250 942.0
קרית אתא 2123.570732 2415.0
ראש פינה 319.470968 226.0

In addition, we calculated the mean MSE presented below. We see that the MSE value is very high and the predictions were not very good.

In [128]:
MSE = mean_squared_error(y_actual, y_pred)
In [129]:
ab_mse = pd.DataFrame(round(MSE, 3), columns=['Mean MSE'], index=['Ada-Boost Regressor'])
ab_mse
Out[129]:
Mean MSE
Ada-Boost Regressor 9177627.127
In [130]:
sum_crimes = crimes.groupby("שם יישוב").sum()
sum_crimes_Q5 = sum_crimes[sum_crimes.index.isin(ab_res.index.values)]

Graphs of the Results

In [131]:
plt = sum_crimes_Q5[sum_crimes_Q5.index.isin(['אילת'])][[2014,2015,2016,2017,2018,2019]].T.plot()
legend(['אילת'[::-1]])
point = pd.DataFrame({'x': [2019], 'y': ab_res[ab_res.index.isin(['אילת'])]['מספר העבירות החזוי']})
plt.vlines(2019, sum_crimes_Q5[sum_crimes_Q5.index.isin(['אילת'])][[2019]], 
         ab_res[ab_res.index.isin(['אילת'])]['מספר העבירות החזוי'], ls='--')
point.plot(x='x', y='y', ax=plt, kind='scatter', label='ערך חזוי'[::-1], color='red')
plt.set_xlabel('Year')
plt.set_ylabel('Number of Crimes')
Out[131]:
Text(0, 0.5, 'Number of Crimes')
In [132]:
plt = sum_crimes_Q5[sum_crimes_Q5.index.isin(['באר שבע'])][[2014,2015,2016,2017,2018,2019]].T.plot()
legend(['באר שבע'[::-1]])
point = pd.DataFrame({'x': [2019], 'y': ab_res[ab_res.index.isin(['באר שבע'])]['מספר העבירות החזוי']})
plt.vlines(2019, sum_crimes_Q5[sum_crimes_Q5.index.isin(['באר שבע'])][[2019]], 
         ab_res[ab_res.index.isin(['באר שבע'])]['מספר העבירות החזוי'], ls='--')
point.plot(x='x', y='y', ax=plt, kind='scatter', label='ערך חזוי'[::-1], color='red')
plt.set_xlabel('Year')
plt.set_ylabel('Number of Crimes')
Out[132]:
Text(0, 0.5, 'Number of Crimes')
In [133]:
plt = sum_crimes_Q5[sum_crimes_Q5.index.isin(['סח\'נין'])][[2014,2015,2016,2017,2018,2019]].T.plot()
legend(['סח\'נין'[::-1]])
point = pd.DataFrame({'x': [2019], 'y': ab_res[ab_res.index.isin(['סח\'נין'])]['מספר העבירות החזוי']})
plt.vlines(2019, sum_crimes_Q5[sum_crimes_Q5.index.isin(['סח\'נין'])][[2019]], 
           ab_res[ab_res.index.isin(['סח\'נין'])][['מספר העבירות החזוי']], ls='--')
point.plot(x='x', y='y', ax=plt, kind='scatter', label='ערך חזוי'[::-1], color='red')
plt.set_xlabel('Year')
plt.set_ylabel('Number of Crimes')
Out[133]:
Text(0, 0.5, 'Number of Crimes')
In [134]:
plt = sum_crimes_Q5[sum_crimes_Q5.index.isin(['קרית אתא'])][[2014,2015,2016,2017,2018,2019]].T.plot()
legend(['קרית אתא'[::-1]])
point = pd.DataFrame({'x': [2019], 'y': ab_res[ab_res.index.isin(['קרית אתא'])]['מספר העבירות החזוי']})
plt.vlines(2019, sum_crimes_Q5[sum_crimes_Q5.index.isin(['קרית אתא'])][[2019]], 
         ab_res[ab_res.index.isin(['קרית אתא'])]['מספר העבירות החזוי'], ls='--')
point.plot(x='x', y='y', ax=plt, kind='scatter', label='ערך חזוי'[::-1], color='red')
plt.set_xlabel('Year')
plt.set_ylabel('Number of Crimes')
Out[134]:
Text(0, 0.5, 'Number of Crimes')
In [135]:
plt = sum_crimes_Q5[sum_crimes_Q5.index.isin(['ראש פינה'])][[2014,2015,2016,2017,2018,2019]].T.plot()
legend(['ראש פינה'[::-1]])
point = pd.DataFrame({'x': [2019], 'y': ab_res[ab_res.index.isin(['ראש פינה'])]['מספר העבירות החזוי']})
plt.vlines(2019, sum_crimes_Q5[sum_crimes_Q5.index.isin(['ראש פינה'])][[2019]], 
         ab_res[ab_res.index.isin(['ראש פינה'])]['מספר העבירות החזוי'], ls='--')
point.plot(x='x', y='y', ax=plt, kind='scatter', label='ערך חזוי'[::-1], color='red')
plt.set_xlabel('Year')
plt.set_ylabel('Number of Crimes')
Out[135]:
Text(0, 0.5, 'Number of Crimes')

6 Police Forces - Recommendation

In this section we attempt to predict the required police forces for each city.

First we create a dataframe that includes only columns with the word עבירות.

In [136]:
cols = [col for col in train_df_scaled.columns if 'עבירות' in col]
offenses = train_df_scaled[cols]

Random Forest Regressor

We loop through the different crime types and using the train data, we predict the number of the specific crime type in 2019.

In [137]:
parameters_grid = {
    'max_depth': [50, 100, 150, 200],
    'min_samples_leaf': [1, 2],
    'min_samples_split': [2, 3],
    'n_estimators': [50, 100, 150],
    'max_features': [5, 10, 15, 20]
}

pred_res = []
MSE = []

# create a loop for each 'עבירות' that is in the dataframe
for i in range(len(offenses.columns)):
    
    # exclude the ith column from the train df
    train_cols = [col for col in train_df_scaled.columns if col not in [offenses.columns[i]]]
    # train df without the ith column
    temp_train = train_df_scaled[train_cols].drop(['סמל יישוב'], axis=1)
    
    # exclude the ith column from the test df
    test_cols = [col for col in test_df_scaled.columns if col not in [offenses.columns[i]]]
    # test df without the ith column
    temp_test = test_df_scaled[test_cols].drop(['סמל יישוב'], axis=1)
    
    # define model
    rf = RandomForestRegressor(random_state=RSEED)
    # define grid search
    grid_search = GridSearchCV(estimator=rf, param_grid=parameters_grid, cv=10, n_jobs=-1)
    # fit the model
    grid_search.fit(temp_train, train_df[offenses.columns[i]])
    # get best model
    best = grid_search.best_estimator_
    # predict using best model
    y_pred = best.predict(temp_test)
    
    # add predictions
    pred_res.append(y_pred)
    
    # add MSE of 
    MSE.append(mean_squared_error(test_df[offenses.columns[i]], y_pred))

Next, we define the average number of work hours required for each crime type.

In [138]:
hours = np.array([20, 5, 15, 10, 35, 0, 25, 40, 30, 1, 1, 2, 0])

Using the function get_overall_hours, we calculate the overall work hours required per city. We calculate this for out train, test and predicted data.

In [139]:
def get_overall_hours(df):
    
    # calculate hours per crime
    i = 0
    for col in offenses.columns.values:
        df[col] = df[col] * hours[i]
        i += 1
    
    # get sum of hours needed
    df['שעות עבודה'] = df[offenses.columns].sum(axis=1)
    
    return df.drop(offenses.columns, axis=1)
In [140]:
train_hours = get_overall_hours(train_df)
test_hours = get_overall_hours(test_df)

rf_res = pd.DataFrame()
rf_res = rf_res.append(pred_res).T
rf_res.columns = offenses.columns

pred_hours = get_overall_hours(rf_res)
In [141]:
res_df = pd.DataFrame([train_hours['שעות עבודה'], test_hours['שעות עבודה'], pred_hours['שעות עבודה']], 
                      index=['train', 'test', 'predicted']).T
In [142]:
res_df.head()
Out[142]:
train test predicted
0 5642.6 5179.0 4166.699421
1 4045.4 4019.0 6955.287804
2 4999.4 3309.0 5718.637155
3 18005.4 20553.0 19159.662016
4 14469.0 13616.0 17432.681259

We divide the cities into 3 clusters in the following matter:

Let $$\text{ratio}_i = \frac{\text{train}_i}{\text{predicted}_i}$$

where $\text{train}_i$ is the average number of hours in the past for city $i$, and $\text{predicted}_i$ is the predicted number of hours in 2019.

If the ratio is smaller than 0.5 (i.e. we predict more hours in 2019), then the police should increase it's forces in city $i$. If the ratio is greater than 0.5 (i.e. we predict less hours in 2019), then the police should decrease it's forces in city $i$. Otherwise (i.e. we predict no significant change), then the police should make no changes in city $i$.

In addition, we divide to data to the same clusters using the test data instead of the predicted data to evaluate our predictions.

In [143]:
for index, row in res_df.iterrows():
    
    # predicted cluster
    if row['train'] / row['predicted'] < 0.5:
        res_df.loc[index, 'pred_cluster'] = 'Increase'
    elif row['train'] / row['predicted'] > 1.5:
        res_df.loc[index, 'pred_cluster'] = 'Decrease'
    else:
        res_df.loc[index, 'pred_cluster'] = 'Keep'
    
    # actual cluster
    if row['train'] / row['test'] < 0.5:
        res_df.loc[index, 'actual_cluster'] = 'Increase'
    elif row['train'] / row['test'] > 1.5:
        res_df.loc[index, 'actual_cluster'] = 'Decrease'
    else:
        res_df.loc[index, 'actual_cluster'] = 'Keep'
        
In [144]:
res_df.head()
Out[144]:
train test predicted pred_cluster actual_cluster
0 5642.6 5179.0 4166.699421 Keep Keep
1 4045.4 4019.0 6955.287804 Keep Keep
2 4999.4 3309.0 5718.637155 Keep Decrease
3 18005.4 20553.0 19159.662016 Keep Keep
4 14469.0 13616.0 17432.681259 Keep Keep
In [145]:
import matplotlib.pyplot as plt

The following plot shows the predicted number of work hours, compared to the true values.

In [146]:
plt.figure(figsize=(8, 6))
ax = plt.gca()
ax.scatter(res_df['test'], res_df['predicted'], color='#b3cde3')
ax.set_yscale('log')
ax.set_xscale('log')
plt.title('Predicted vs. True Number of Hours')
plt.xlabel('Actual')
plt.ylabel('Predicted')
plt.grid(alpha=0.5, color='lightgrey')
plt.show()

The following confusion matrix presents the predicted classification of the cities, compared to the true classification.

In [147]:
matrix = confusion_matrix(res_df['actual_cluster'], res_df['pred_cluster'])

plt.figure(figsize=(8, 6))
sns.heatmap(matrix, annot=True, cmap='Blues',
            xticklabels=['Increase', 'Decrease', 'Keep'], yticklabels=['Increase', 'Decrease', 'Keep'])
plt.xlabel("True", fontsize=12)
plt.ylabel("Predicted", fontsize=12)
plt.title("Confusion Matrix", fontsize=14)
plt.show()
In [148]:
acc = pd.DataFrame([round(accuracy_score(res_df['actual_cluster'], res_df['pred_cluster']), 3)],
                  index=['Classification'], columns=['Accuracy'])
acc
Out[148]:
Accuracy
Classification 0.894

We see that most cities were correctly classified (accuracy score 0.894). However, note that the data is quite imbalanced with much more observations in the 'Keep' class.


7 Random Forset with Clusters

In this section, we run the Random Forest algorithm from section 4 again. This time, we add the clusters from section 3 (K-Means with $K=3$).

First, we add the clusters to the train and test datasets.

In [149]:
train_df_scaled['Cluster'] = kmeans_df3['Cluster']
test_df_scaled['Cluster'] = kmeans_df3['Cluster']

Then, we create a dataframe that includes only columns with the word עבירות.

In [150]:
cols = [col for col in train_df_scaled.columns if 'עבירות' in col]
offenses = train_df_scaled[cols]
offenses.head()
Out[150]:
עבירות בטחון עבירות כלכליות עבירות כלפי המוסר עבירות כלפי הרכוש עבירות מין עבירות מנהליות עבירות מרמה עבירות נגד אדם עבירות נגד גוף עבירות סדר ציבורי עבירות רשוי עבירות תנועה שאר עבירות
0 -0.327640 -0.403828 -0.566207 -0.409107 -0.404847 -0.280400 -0.386804 -0.197029 -0.267263 -0.343166 -0.391465 -0.320557 0.091901
1 -0.236044 -0.437452 -0.414116 -0.429182 -0.520959 -0.280400 -0.430506 -0.625057 -0.397166 -0.427081 -0.344954 -0.396237 -0.321652
2 -0.678757 -0.437452 -0.536634 -0.239811 -0.357347 0.006521 -0.377093 -0.563910 -0.454725 -0.467093 -0.530997 -0.497143 -0.321652
3 3.824702 0.235029 -0.023327 -0.060375 -0.367902 0.006521 -0.097081 1.453936 0.050630 0.298419 0.073642 0.814642 0.091901
4 -0.350539 -0.202084 0.013640 -0.147713 0.133492 -0.136939 -0.207144 -0.013588 0.005745 -0.002229 0.748046 0.259656 0.919005

Next, we select the required rows from the test dataset.

In [151]:
cities_sign = ([1139, 7400, 507, 2640, 43, 3780])
rf_test = test_df[test_df['סמל יישוב'].isin(cities_sign)]
rf_test
Out[151]:
סמל יישוב סך הכל אוכלוסייה 2018 ערבים שנת ייסוד קואורדינטות גובה דת יישוב_2.0 דת יישוב_3.0 דת יישוב_4.0 דת יישוב_אחר ... עבירות מין עבירות מנהליות עבירות מרמה עבירות נגד אדם עבירות נגד גוף עבירות סדר ציבורי עבירות רשוי עבירות תנועה שאר עבירות שעות עבודה
32 3780.0 56746.0 2.0 1985 2.107362e+09 568.0 0 0 0 0 ... 1680.0 0.0 700.0 40.0 3990.0 230.0 3.0 0.0 0.0 11093.0
94 507.0 10029.0 10008.0 1968 2.155976e+09 34.0 1 0 0 0 ... 280.0 0.0 400.0 120.0 2160.0 183.0 2.0 2.0 0.0 6492.0
104 1139.0 46124.0 1507.0 1964 2.276276e+09 69.0 0 0 0 0 ... 1645.0 0.0 2400.0 120.0 8550.0 563.0 15.0 18.0 0.0 21536.0
118 43.0 1599.0 65.0 1896 2.541580e+09 260.0 0 0 0 0 ... 35.0 0.0 50.0 0.0 150.0 9.0 0.0 2.0 0.0 596.0
135 7400.0 217244.0 535.0 1929 1.872569e+09 -2.0 0 0 0 0 ... 3780.0 0.0 9500.0 240.0 40980.0 2354.0 18.0 52.0 0.0 101339.0
177 2640.0 56344.0 31.0 1950 1.961967e+09 29.0 0 0 0 0 ... 1960.0 0.0 1925.0 80.0 8820.0 572.0 6.0 8.0 0.0 21661.0

6 rows × 217 columns

Random Forest Regressor

We loop through the different crime types and using the train data, we predict the number of the specific crime type in 2019.

In [152]:
parameters_grid = {
    'max_depth': [50, 100, 150, 200],
    'min_samples_leaf': [1, 2],
    'min_samples_split': [2, 3],
    'n_estimators': [50, 100, 150],
    'max_features': [5, 10, 15, 20]
}

pred_res = []
MSE = []

# create a loop for each 'עבירות' that is in the dataframe
for i in range(len(offenses.columns)):
    
    # exclude the ith column from the train df
    train_cols = [col for col in train_df_scaled.columns if col not in [offenses.columns[i]]]
    # train df without the ith column
    temp_train = train_df_scaled[train_cols].drop(['סמל יישוב'], axis=1)
    
    # exclude the ith column from the test df
    test_cols = [col for col in test_df_scaled.columns if col not in [offenses.columns[i]]]
    # test df without the ith column
    temp_test = test_df_scaled[test_cols].drop(['סמל יישוב'], axis=1)
    
    # define model
    rf = RandomForestRegressor(random_state=RSEED)
    # define grid search
    grid_search = GridSearchCV(estimator=rf, param_grid=parameters_grid, cv=10, n_jobs=-1)
    # fit the model
    grid_search.fit(temp_train, train_df[offenses.columns[i]])
    # get best model
    best = grid_search.best_estimator_
    # predict using best model
    y_pred = best.predict(temp_test)
    
    # add predictions for the required cities
    pred_res.append(y_pred[rf_test.index])
    
    # add MSE of 
    MSE.append(mean_squared_error(test_df[offenses.columns[i]], y_pred))

The following table presents the predicted crime numbers for each test city.

In [153]:
rf_res = pd.DataFrame()
rf_res = rf_res.append(pred_res).T
rf_res.columns = offenses.columns
rf_res.index = ['כרמיאל','נתניה','כפר יאסיף','ראש העין','מטולה','ביתר עילית']
rf_res
Out[153]:
עבירות בטחון עבירות כלכליות עבירות כלפי המוסר עבירות כלפי הרכוש עבירות מין עבירות מנהליות עבירות מרמה עבירות נגד אדם עבירות נגד גוף עבירות סדר ציבורי עבירות רשוי עבירות תנועה שאר עבירות
כרמיאל 478.52 15.559000 1779.278500 3325.306667 1122.38 0.0 1046.900000 111.600000 5872.64 377.285543 3.655867 8.660 0.0
נתניה 336.08 6.681667 880.979190 1508.840000 131.60 0.0 415.500000 49.860698 2805.96 206.792556 2.464533 6.416 0.0
כפר יאסיף 466.48 36.192333 3858.416000 6472.413333 1617.28 0.0 2162.466667 122.969397 10459.88 697.878093 10.269867 14.116 0.0
ראש העין 46.56 1.918667 196.328667 487.226667 65.24 0.0 91.866667 11.795556 474.72 35.620711 0.569600 1.916 0.0
מטולה 606.36 85.084667 7874.826889 31808.386667 4125.10 0.0 9855.933333 349.841651 36979.76 2388.756911 21.303667 43.264 0.0
ביתר עילית 202.84 19.610000 2781.766381 7911.160000 1155.84 0.0 2228.133333 99.072508 8613.20 648.133644 7.333400 16.316 0.0

Now, we select the second most frequent crime according to the predictions and compare it to the actual data of 2019.

In [154]:
pred_second = rf_res.apply(lambda row: row.nlargest(2).values[-1], axis=1)
pred_ind = rf_res.apply(lambda row: row.nlargest(2).idxmin(), axis=1)
In [155]:
actual_second = rf_test[cols].apply(lambda row: row.nlargest(2).values[-1], axis=1)
actual_ind = rf_test[cols].apply(lambda row: row.nlargest(2).idxmin(), axis=1)
In [156]:
res_df = pd.DataFrame(pred_ind, columns=['עבירה חזויה'])
res_df['עבירה אמיתית'] = actual_ind.values
res_df['מספר העבירות החזוי'] = pred_second
res_df['מספר העבירות האמיתי'] = actual_second.values
In [157]:
res_df
Out[157]:
עבירה חזויה עבירה אמיתית מספר העבירות החזוי מספר העבירות האמיתי
כרמיאל עבירות כלפי הרכוש עבירות כלפי הרכוש 3325.306667 2480.0
נתניה עבירות כלפי הרכוש עבירות כלפי המוסר 1508.840000 1815.0
כפר יאסיף עבירות כלפי הרכוש עבירות כלפי הרכוש 6472.413333 4740.0
ראש העין עבירות נגד גוף עבירות נגד גוף 474.720000 150.0
מטולה עבירות כלפי הרכוש עבירות כלפי הרכוש 31808.386667 33790.0
ביתר עילית עבירות כלפי הרכוש עבירות כלפי הרכוש 7911.160000 6140.0

We see that 5 out of the 6 predictions were accurate, which is a good result.

In addition, we calculated the mean MSE presented below.

In [158]:
rf_mse = pd.DataFrame(round(np.mean(MSE), 3), columns=['Mean MSE'], index=['Random Forest Regressor'])
rf_mse
Out[158]:
Mean MSE
Random Forest Regressor 712041.782

The MSE without clusters was much smaller which is better.